***DRAW FIGURE 4, Caloric Potential Cell***

version 15.1

clear
use "worms_replication.dta"
set more off

*keep only covered dioceses*
drop if covered == 0




egen year_cal = group(year calories)
collapse (mean) religious_bishop year calories, by(year_cal)

rename religious_bishop avg_religious_c




*calculate and record the 40 year moving average share of religious bishops for high and low caloric potential dioceses*
gen moving_avg_c = .
quietly forval i = 800/1289{
	sum avg_religious_c if year >= `i'-20 & year <= `i'+20 & calories == 1 
	replace moving_avg_c = r(mean) if year == `i' & calories == 1
	sum avg_religious_c if year >= `i'-20 & year <= `i'+20 & calories == 0
	replace moving_avg_c = r(mean) if year == `i' & calories == 0
	}

label var avg_religious_c "Share Religious Bishops"	
label var moving_avg_c "Share Religious Bishops"
label var year "Year"

*For each year and each of high and low caloric potential plot both the share of religious bishops (avg_religious_c) and the moving average*

twoway (scatter avg_religious_c year if  calories == 0  & year >=800 & year <=1309, mcolor(black) lcolor(black))(line moving_avg_c year if calories == 0  & year >=800 & year <=1309, mcolor(black) lcolor(black))(scatter avg_religious_c year if   calories ==1  & year >=800 & year <=1309, mcolor(gs12) lcolor(gs12))(line moving_avg_c year if calories ==1  & year >800 & year <=1309, mcolor(gs12) lcolor(gs12) xline(1122, lpattern(dash) lcolor(black)) xline(1309, lpattern(dash) lcolor(black))), ///
   graphregion(color(white))  bgcolor(white) scheme(lean1) xlab(800 900 1000 1100 1200 1300 )  ylab(, nogrid) legend(order(1 "poor diocese" 2 "poor diocese" 3 "wealthy diocese" 4 "wealthy diocese")) title("Caloric Potential") title("Caloric Potential")

